#if !defined(_Qhull_Link_)
#define _Qhull_Link_

#include <src/data/all.h>
#include <src/link/qhull/qhull_a.h>
#include <src/math/Operators/ParticlesOperators.h>

namespace QHULL 
{
	typedef particle<OS_double>										p_array_1;
	typedef particles_list<p_array_1>								p_array_2;

	typedef p_array_2												particles_list;
	typedef particles_list::Base									Base;
	typedef Base::value_type										particle_type;
	typedef Base::iterator											comp_iterator;
	typedef particle_type::iterator									part_iterator;
	typedef particles_list::indice_type								indice_type;
	typedef particles_list::indice_type								indice_type;
	typedef graph_base<particles_list>								graph_base_type;
	typedef graph_base_type::weight_set_type						weight_set_type;
	typedef graph_base_type::graph_type								edges_set_type;
	typedef graph_base_type::edge_type								edge_type;

	template<class A, class R, std::size_t Dim = A::Dim>
	class qh_edges_set_type 
		: public ParticlesOperators<A,R>
	{
		OS_STATIC_CHECK("false"); // this class should never be compiled. Specialization to 1D and ND below.
	};



	template<class A, class R>
	class qh_edges_set_type<A,R,1> 
		: public ParticlesOperators<A,R>
	{
	public:
		typedef  Particle_Operator										Family_Operator;
		typedef  ParticlesOperators										Base;
		typedef  typename Base::Base_Cast								Base_Cast;
		typedef qh_edges_set_type										type;

		qh_edges_set_type( ) {};
		qh_edges_set_type(const Base_Cast &first) : Base(first)
		{};
		qh_edges_set_type(const qh_edges_set_type & from) :  Base(from) {
		};
		qh_edges_set_type & operator = (const qh_edges_set_type & from) {
			Base::operator =(from);
			return *this;
		};
		virtual void update() {
			std::sort(get_state_variable()->begin(), get_state_variable()->end());
		};
		virtual ~qh_edges_set_type(){};
		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,weight_set_type>(&weight_set_);
		};

	};

	template<class A, class R>
	class qh_edges_set_type<A,R,2> 
		: public ParticlesOperators<A,R>
	{
	public :
		typedef  Particle_Operator										Family_Operator;
		typedef  ParticlesOperators										Base;
		typedef  typename Base::Base_Cast								Base_Cast;
		typedef qh_edges_set_type										type;

		qh_edges_set_type( ) {};
		qh_edges_set_type(const Base_Cast &first) : Base(first)
		{};
		qh_edges_set_type(const qh_edges_set_type & from) :  Base(from) {
		};
		qh_edges_set_type & operator = (const qh_edges_set_type & from) {
			Base::operator =(from);
			return *this;
		};
		virtual ~qh_edges_set_type(){};
		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,weight_set_type>(&weight_set_);
		};
		virtual void update() {
			construct<weight_set_type>();
		};
		template <class T>
		void construct()
		{
			OS_STATIC_CHECK(false); // you must implement the Delaunay graph for your own weight set type
		};
		template <>
		void construct<weight_set_default_def<base_value_type>>() 
		{
			weight_set_.clear();
			OS_int numpoints = get_state_variable()->size(2), dim = get_state_variable()->size(1);
			weight_set_.resize(numpoints,numpoints,false);
			if (numpoints == 0)
				return;
			OS_DYNAMIC_CHECK(dim >1, "Non correct dimension in Delaunay computations");
			coordT* points = new coordT[numpoints*dim]; /* array of coordinates for each point */
			{
			  OS_int j,k;
			  coordT *point;

			  for (j=0; j<numpoints; j++) {
				point= points + j*dim;
				for (k= 0; k < dim; k++) {
				  double breakpoint = (*get_state_variable())[j][k];
				  point[k]= breakpoint;
				}
			  }
			}

			boolT ismalloc= False;    /* True if qhull should free points in qh_freeqhull() or reallocation */
			FILE* outfile = NULL;    
		/* output S qh_produce_output()
					   use NULL to skip qh_produce_output() 
					   use stdout for verbose 
		*/
			FILE* errfile = stderr;    
									/* error messages S qhull code. */

									  /*qhull s Tcv : parametres de qhull. Voir doc"*/
									  /*option Pp : pour desactiver le check. Voir doc"*/
									  /* "qhull v Qbb Qv Qx" marche mais est lent*/

	/* duplicated in qdelau_f.htm and qdelaun.htm */


		  OS_int exitcode= qh_new_qhull 
			  (dim, 
			  numpoints, 
			  points, 
			  ismalloc, 
			  "qhull d v Qbb", 
			  outfile, 
			  errfile);

		  OS_DYNAMIC_CHECK (exitcode == 0, "Delaunay Algorithm (qhull) failed" );
#if defined(_DEBUG) 
		  qh_check_points ();
#endif 

			qh_vertexneighbors();
			vertexT *vertex, **vertexp;
			facetT *neighbor,**neighborp;
			qh vertex_visit++;
			FORALLvertices
			{
	//			std::cout << std::endl << " ------------visiting vertex : " << vertex->visitid << "--point id " << qh_pointid (vertex->point) << std::endl; 
			  vertex->visitid = qh vertex_visit;
			  unsigned OS_int idA = qh_pointid (vertex->point);
			  FOREACHneighbor_(vertex)
			  {
				if (neighbor->upperdelaunay)
				{
	//				std::cout << "     --- Delaunay neighbors : " << neighbor->visitid << std::endl; 
				}
				else
				{
	//				std::cout << "     --- visiting vertex neighbors : " << neighbor->visitid << std::endl; 
				  neighbor->visitid = qh vertex_visit;

				  vertexT *localvertex,**localvertexp;

				  FOREACHsetelement_(vertexT,neighbor->vertices,localvertex)
				  {
					unsigned OS_int idB = qh_pointid (localvertex->point);
					if(localvertex->visitid != qh vertex_visit)
					{
	//						std::cout << "        " << idA << "----" << idB << " ... inserted S vertex " << vertex->visitid; 
	//						std::cout << std::endl;
						unsigned OS_int idB = qh_pointid (localvertex->point);

						weight_set_(idA,idB) = 1.;
						weight_set_(idB,idA) = 1.;
					}
				  }
				}
			  }
			}

		  qh_freeqhull( true );
		  OS_int curlong, totlong;
		  qh_memfreeshort (&curlong, &totlong);

		  if (curlong || totlong)
			fprintf (qh ferr, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
		  delete [] points;
		}
//////////////////////////////////////////
	};
}; // namespace QHUll
	
////////////////////////////////////////

#endif